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We calculate pion vector and scalar form factors in two-flavor lattice QCD and study the 
chiral behavior of the vector and scalar radii {r'^)v,s- Numerical simulations are carried out 
on a 16'^ X 32 lattice at a lattice spacing of 0.12 fm with quark masses down to ~ TOs/6, 
where rrig is the physical strange quark mass. Chiral symmetry, which is essential for a direct 
comparison with chiral perturbation theory (ChPT), is exactly preserved in our calculation 
at finite lattice spacing by employing the overlap quark action. We utilize the so-called all-to- 
all quark propagator in order to calculate the scalar form factor including the contributions 
of disconnected diagrams and to improve statistical accuracy of the form factors. A detailed 
comparison with ChPT reveals that the next-to-next-to-leading-order contributions to the 
radii are essential to describe their chiral behavior in the region of quark mass from TOs/6 to 
ms/2. Chiral extrapolation based on two-loop ChPT yields (r'^) y = 0.409(23)(37) fm^ and 
{i^^)s = 0.617(79)(66) fm^, which are consistent with phenomenological analysis. We also 
present our estimates of relevant low-energy constants. 
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I. INTRODUCTION 



Recent algorithmic improvements allow us to perform large-scale simulations of unquenched 
QCD in the chiral regime with various lattice discretizations. Calculation of phenomenologically 
important quantities has then become feasible. In order to make reliable prediction for physical 
observables, it is crucial to examine the consistency of their chiral behavior with expectations 
from chiral perturbation theory (ChPT). Lattice QCD with exact chiral symmetry is the cleanest 
framework for this purpose, while conventional lattice actions may distort chiral behavior of 
observables by their explicit symmetry breaking, especially when one goes beyond the next- 
to-leading order (NLO) in ChPT. The JLQCD and TWQCD collaborations embarked on the 



simulations with exact chiral symmetry 



QCD [5|]. For other physics results from this project, see [a, lZ|, 
The pion vector form factor Fy(g^) defined by 

{7t{p')\VMp)) = [p + p')Fv{q^), 



4] employing the overlap quark action 

0,B- So 

cc 

i 



far, 



we have performed a detailed study of the pion mass and decay constant in two-flavor lattice 



101, 



11 



{v-p')\ 



(1) 



provides a simple testing ground for the consistency between lattice calculation and ChPT. 
This is one of the fundamental quantities to characterize the low-energy dynamics of pions: for 
instance, it is related to the charge radius of pion 

dFv{q^)\ 

Ig2=0 



6 



dq^ 



(2) 



able up to two loops, namely to the next-to- 
, Q, Q and A^/ = 3 



14 



The chiral expansion of Fv{q'^) and (r^)y is avai 
next-to-leading-order (NNLO), both for Nj = 2 

experimental data based on two- loop ChPT have led to precise estimates of (r^)y [ifi], [l^, which 
can also be used as a benchmark of lattice calculations 



171 . llSl l. Analyses of 
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From previous lattice studies 
periment has not 
experiment 24. [26 



19|, 



21 



22 



24 



251 . I26l . |27I |. the consistency with ex- 
agreement with 



3een established convincingly: some of them reported gooc 
. [27I. whereas others underestimated (r^)y significantly 



201, 



22 



23 



251 . This 



is possibly due to systematics in the parametrization of the q^ dependence of Fv{q^) and to 
the chiral extrapolation of {r'^)v- For instance, the NNLO chiral corrections are fully taken 
into account in the ChPT analyses but not in most of the previous lattice studies. They could 
significantly modify the chiral behavior of (r-^)y at up and down quark masses larger than their 
physical value, as demonstrated in our report [28] and more recently in Ref. with a different 
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lattice discretization. Our simulations with exact chiral symmetry enable us to perform more 
stringent test of the chiral behavior of lattice data using two-loop ChPT without suffering from 
the distortion due to explicit chiral symmetry breaking present in other frameworks. 

The chiral behavior of the scalar form factor Fs^q"^) defined by 

{Ap')\S\7:{p)) = Fs{q') (3) 
is another interesting subject, since i) its radius 



/ 2v a 9Fsiq^) 



(4) 

<?2=0 



dq 

provides a determination of the low energy constant (LEG) I4 alternative to that with the 
decay constant -Ftt, and ii) {r'^)s has 6 times enhanced chiral logarithm compared to {r'^)v B^nd 
thus may offer an opportunity to clearly identify the one-loop chiral logarithm. Since there 
are no experimental processes directly related to Fsiq"^), its direct determination is possible 
only through lattice QCD. It is however difficult to evaluate disconnected correlation functions 
with conventional point-to-all quark propagator, which flows from a fixed lattice site to any site. 
There have been only a few calculations of Fs{q'^) 23, l25(], and the contributions of disconnected 
diagrams were ignored in these studies. 

In this paper, we calculate the pion vector and scalar form factors Fy^siq^) in two- flavor 
QCD and study the chiral behavior of the radii {r'^)v,s- For a detailed comparison with two- 
loop ChPT, we preserve chiral symmetry by employing the overlap quark action and simulate 
up and down quarks with masses as low as m~ms/6. The scalar form factor Fs{q'^) is evaluated 
including the contribution of disconnected diagrams by using the all-to-all quark propagator, 
which contains propagations from any lattice site to any site. The all-to-all propagator is also 
helpful to substantially improve statistical accuracy of Fv^siq"^)- O ur p reliminary analyses based 
on one- and two-loop ChPT have been reported in Refs. [23] and [28l |, respectively. 

This paper is organized as follows. We introduce our simulation method in Sec. HH Calcula- 
tion of the form factors from pion correlators is presented in Sec lllli We parametrize their q'^ 
dependence in Sec. llVi Section |V] is devoted to a detailed description of our chiral extrapolation 
of the radii. Finally, our concluding remarks are given in Sec. IVII 
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II. SIMULATION METHOD 



A. Configuration generation 



We calculate pion form factors in QCD with dynamical up and down quarks with a de gen erate 



^ and 



mass parameter. Numerical simulations are carried out with the Iwasaki gauge action 
the overlap quark action [s, 4|, which has exact chiral symmetry at finite lattice spacings 3l| 
Its Dirac operator is given by 



D{m) 



m 



mo + — ] + [mo 



m 



75 sgn [Hy/i-mo]] 



(5) 



where m is the quark mass and H^j^{—mQ)=^^D^j^{—mo) is the Hermitian Wilson-Dirac opera- 



tor. We set the mass parameter of this kernel 
of the overlap Dirac operator is confirmed [32 



operator to —mo = —1.6, with which the locality 



33l |. Because of the sign function sgn[//w(— ""^o)] 



in D{m), the overlap action is discontinuous when H^{—mo) develops zero eigenvalue (s). The 
commoiily-used Hybrid Monte Carlo (HMC) algorithm can be modified to deal with this discon- 



tinuity [3j] but turned out to be very costly. In order to carry out high-statistics simulations 
we suppress (near-)zero modes of H^^{—mo) by introducing an auxiliary determinant 

dei[H^{-mof] 



A 



w 



(6) 



into the Boltzmann weight 



35], 



36( 1 . We note that this can be considered as an O(a^) modification 
of the gauge action, and hence does not change the continuum limit of the theory. The parameter 
/i is tuned to 0.2 so as to minimize lattice artifacts induced by Aw 32, l36|]. We refer readers to 
Ref. [33] for further details of our simulation method. 

An important property of the determinant Aw is that it fixes the global topological charge Q 
of the gauge field during continuous updating of the gauge configuration in the HMC algorithm. 
Note, however, that local topological fluctuations are present, and the topological susceptibility 
calculated in Ref. shows expected behavior as a function of sea quark mass. As shown in 



Refs. 



37 



38( 1 ■ the effect of fixed global topology can be considered as a finite volume effect 



and is suppressed by the inverse of the space-time volume V . Furthermore the effect can be 
systematically corrected by investigating Q dependence of physical observables of interest 



Our gauge ensembles are generated on a 16'^ x 32 lattice. The bare gauge coupling is 13 = 
Q/g^ = 2.30 where the lattice spacing determined from the Sommer scale rg = 0.49 fm [391] is 
a = 0.1184(21) fm. In the trivial topological sector (5 = 0, we take four values of bare up and 
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down quark masses, m = 0.015, 0.025, 0.035 and 0.050, which cover a range [m^/G, ms/2]. At 
each m, we calculate pion form factors using 100 independent configurations separated by 100 
HMC trajectories. In order to study the effect of fixed topology, we also simulate non-trivial 
topological sectors Q = —2 and —4 at m = 0.050 with statistics of 50 independent configurations. 



B. Construction of all-to-all quark propagator 

The conventional method to calculate hadron correlators employs the so-called point-to-all 
quark propagator which flows from a fixed lattice site to any site. This is however not suitable to 
calculate disconnected diagrams which involve quark loops starting from and ending at arbitrary 
lattice sites. In this work, therefore, we construct all-to-all quark propagator that contains the 
quark propagating from any lattice site to any site along the strategy proposed in Ref. [4^. 

Let us consider a decomposition of the quark propagator using eigenmodes of the overlap 
operator 

D-\x,y) = 5^-Ln«(x)nWt(y), (7) 

k 

where A^'^) and u^^\x) represent k-th lowest eigenvalue and eigenvector oi D, respectively. Note 
that the eigenvalues are ordered by their absolute values. Color and spinor indices are suppressed 
for simplicity. It is expected that low-lying modes dominate low-energy dynamics of pions 
including their form factors. We evaluate these low-mode contributions to the propagator exactly 
as 

(I)-i)io„(x,y) = ^ n('=)(x)n('=)t(y) (iVe = 100) (8) 
k=i 

using 100 low-lying modes for each gauge configuration. Note that the overlap operator is normal 
and we do not have to distinguish left and right eigenvectors. 

The contribution of higher modes is estimated stochastically by the noise method with the 
dilution technique 40]. We prepare a single Z2 noise vector rj^'^^x) for each configuration, 
and dilute it into = 3 x 4 x A't/2 vectors, which have non-zero elements only for a single 
combination of color and spinor indices and at two consecutive time-slices. The high-mode 
contribution can be estimated as 



(I?-')high(x,y) = ^x('^)(x)77Wt(y) (9) 

d=l 



6 



by solving a linear equation for each diluted source 

^ y) x^"^ (y) = Y^i^^y " ^iow(x, y)) t?^'^) (y) (d = 1, iV^), (10) 
y y 

where d is an index to represent the dilution and 

flow(x,y) = Y.u^^\x)u^^)\y) (11) 
k=l 

is the projector to the eigenspace spanned by the low-modes. 

By combining Eqs. dH]) and ([9]), the all-to-all quark propagator can be expressed by a matrix 

D-\x,y) = ^?;('=)(x)u;(*^)t(y) (12) 

k=l 

constructed from the following two sets of vectors v and w 

{«;«,..., = {n«,...,n(^=),r?«,...,,?(^'')}, (14) 



where = Ne + Na- 

C. Measurement of pion correlators 

Using the all-to-all propagator (jl2p . the pion two-point function with a temporal separation 
At and a spatial momentum p can be expressed as 

, Nt 

* t=l k,l=l 

Here C'r,</>(i; p) is constructed from the v and if vectors as 

4'i\f^P) = j;<A(r)u;WUx + r,t)rt;«(x,t)e-P-, (16) 

x,r 

and represents the meson field with a Dirac matrix T and a momentum p at a temporal coordi- 
nate t. We use a local i;^'i(r) = (5r,o and an exponential form 0s(r) = exp[— 0.4|r|] for the smearing 
function 0(r) in this study. Note that the source point x is averaged over spatial volume. 

Pion form factors are extracted from three-point functions shown in Fig. [H which can also 
be calculated from the meson fields Or,(f,{t, p) as 

Nt N, 

Nt 



, Nt N, 

C^T^ ( At , At'; p, p') = ^ Y + + At'; p') O^^;^] (t + At; p - p') 



t=l k,l,m=l 





<r> (q^=0) 




(conn) (disc) (vev) 

FIG. 1: Connected (left-most diagram) and disconnected three point functions (middle diagram). Note 
that Fs{0) receives a contribution from the right-most diagram due to the non-zero vacuum expectation 
value (VEV) of the scalar operator S. 



, Aft 

* t=l k,l=l 

k=l 



, Aft 

t 4- 1 7, ] -1 

' Aft Ar„ 



i=l k,l=l 

Nt Af„ 




X > 

\ Nt 

\ ^ t'=lk=l 

where (• • • )conf represents the ensemble average. We denote the temporal separation and spatial 
momentum for the initial (final) meson by At and p (At' and p'), respectively. 

We prepare the v and w vectors on the IBM Blue Gene/L at KEK, which has the peak speed of 
57.3 TFLOPS. We employ the implicitly restarted Lanczos algorithm to calculate low-modes of 
D. The computational cost of this step is roughly 0.6 TFLOPS -hours per configuration. Solving 
Eq. (llOp is the most time-consuming part in our measurement, since it requires Nt/2 times 
more inversions than the conventional measurement of two-point functions with the point-to-all 
propagator. We use the conjugate gradient (CG) algorithm accelerated by the relaxed stopping 



condition 4l[| and by the preconditioning using the 100 low-modes. The resulting CPU cost 
of our CG solver is ~ 1.7 TFLOPS • hours/conf. The calculation of the meson field Op'^^'^(t,p) 
needs much less CPU time than the above two steps: it is about 0.2 GFLOPS • hours/conf for a 
single choice of (m,p,r, (p). Once 0{.^;'^(t, p) is prepared, we can calculate all of the connected 
and disconnected pion correlators with small additional cost. These calculations are carried out 
on the Hitachi SRllOOO with the peak speed of 2.15 TFLOPS and workstations at KEK. 

It is advantageous that we do not have to repeat the time consuming preparation of v and 
w vectors to calculate meson fields C'r,0(t,p) with different choices of p, F and 0. In order to 



8 



^ 



m = 0.025, lpl = sqrt(2), lp'l=0, Ar = 7. At' = 1 
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FIG. 2: Statistical fluctuation of connected three point functions 



(conn) 

V^TT .(f) s<f) s 



(Ai, Ai';p,p') (left panel) 



and C, 



(conn) 
7rl7r,0s( 



(Ai, Ai';p, p') (right panel) with At = At' = 7, |p| \/2 and |p'| = 0. Open triangles 
are jackknife data with a single choice of the source location (x, t) and momentum configuration (p,p'), 
whereas the filled squares are averaged over (x, t) and (p, p') corresponding to the same value of . Each 
data is normalized by its statistical average. 



simulate various values of g^, we take 27 choices of the spatial momentum p with |p| < \/3 for 
the initial and final pions, and 33 choices with |q| <2 for the momentum transfer q. Note that 
the spatial momentum is shown in units of 27r/L in this article. This setup covers the region of 
momentum transfer — 1.7< [GeV^] <0. 

The use of the all-to-all propagator also improves the statistical accuracy of the pion corre- 
lators by averaging over the locations of source operator*, namely (x, t) in Eqs. (|16p - (jl9p . as 
well as over momentum configurations (p, p') corresponding to the same value of q^. Figure [2] 
compares the statistical fluctuation of connected pion correlators with a certain choice of At, 
At' and (f'. We observe that averaging over the source locations and momentum configurations 
remarkably reduces the statistical error of the correlators and leads to very accurate results for 
the form factors. 

* We note in passing that only the low-mode contribution is averaged over the locations of the source operator in 
the so-called low-mode averaging (LMA) method [13, 5^. The high-mode contribution in LMA is estimated for 
a fixed source location using the point-to-all propagator in contrast to our method with the all-to-all propagator. 
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III. DETERMINATION OF PION FORM FACTORS 
A. Ratio method for vector form factor 



In the limit of large temporal separations among pion operators and vector current, two- and 
three-point functions are dominated by the contribution from the ground state 

^(conn).. . / ,s -^TT.t^'dp'D* ■^7r,(^(|p|) 1 , , /mt^i / n\ 

C^,^ (At,At;p,p) ^^^^^ 4i^.(|p1)i^.(|p|) ^(-(p^\^Mp)) 

xe--^-(|p'l)^*'e~'^'^(|p')^*, (21) 

where Z7r,(^(|p|) = {'^ip)\0^r^^(i,) is the overlap of the interpolating field 075,,/) to the physical 
state, and Zy is the renormalization factor for the vector current. For a precise determination 
of the matrix element (7r(p')|V^|7r(p)), it is advantageous to take a ratio of appropriately chosen 
correlators in order to cancel out the exponential damping factors e^^'^'-'P^ ''^^*^'' and other 



unnecessary factors Z^^^(,) and Zy [4J|. In this study, we use the following ratio to calculate an 

effective value of the vector form factor Fv{q'^) 

p .A, A,/.„2^ _ 2M^ Rv{At,At';\p\,\p'\,q^) 

''^ ' ^ " ii^.(|p|)+i?.(|p'|) i?y(At, At'; 0,0,0) ' 

„ /A+ A+/ I I I /| 2^ 1 C'i-74,r ■'(At, At'; p, p') 
i?y(At,At;|p|,|p |,g ) = 2^ ^^ (/\ff-w\ ' 

IpI'Ip'I fixed^,|p'| '-^'^'^.'A.^iV^*' P) '-vr7r,</,,0,(At , p j 

where (l/-^|p|,|p'|) X^fixed |p| |p'| represents the average over momentum configurations corre- 
sponding to the same value of q'^. 

The kinematical factor in Eq. (|22p involves energies of the initial and final pion states. Ef- 
fective values of i?7r(|p|) at our largest and smallest quark mass are plotted in Figs. [3] and [H 
Our data show clear signals up to the largest spatial momentum |p| < \/3, since their statistical 
errors are greatly reduced by the use of the all-to-all propagators. We determine i?7r(|p|) by 
a single-cosh fit to (7,^71, 0^03 (At; p) with a fit range At G [At^[^, Nt/2]. The minimum tem- 
poral separation Atmin is chosen by inspecting Atmin dependence of the fit result. Numerical 
results for i?7r(|p|) are summarized in Table [H They are consistent with the dispersion relation 



-^7r(|p|) = \/ + p2, which is commonly assumed in previous studies to estimate ii^T^dpl) for 
IpI > 0. In this work, however, we use the measured value of -Ejrdpl) in order not to underestimate 
uncertainty in £'^(|p|) and hence Fy (At, At'; g^). 

We extract the vector form factor Fv{q^) by a constant fit to the effective value 
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FIG. 3: Effective value of pion energy E'.n.dpl) at (Q, to) = (0, 0.050). Circles, triangles and squares show 
results from CTrTr,^,^, (At; p), C^^r, (At; p) and CrTr^^.^, (At; p), respectively. 
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FIG. 4: Effective value of pion energy i?7r(|p|) at (Q, m) = (0, 0.015). 



TABLE I: Fit results for pion energy i?7r(|p|)- 



Q TO 


IphO 


|P| = 1 


IpHx/2 


|phV3 





0.015 


0.1722(20) 


0.4236(88) 


0.554(27) 


0.710(17) 





0.025 


0.2193(16) 


0.4499(58) 


0.5791(96) 


0.694(16) 





0.035 


0.2610(15) 


0.4706(33) 


0.6091(59) 


0.721(11) 





0.050 


0.3128(12) 


0.5005(21) 


0.6292(28) 


0.732(11) 


-2 


0.050 


0.3124(15) 


0.4972(34) 


0.6257(67) 


0.738(14) 


-4 


0.050 


0.3155(17) 


0.4994(29) 


0.6353(51) 


0.705(10) 



Fv{At,At';q'^). Examples of Fy (At, At'; g^) are plotted in Figs. [5] and H The all-to-all quark 
propagator enables us to change At and At' independently in contrast to previous lattice studies 
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with At 



Effective value of vector form factor Fy (At, At'; q^) at (Q, m) = (0, 0.050). We only plot data 
h At' < Nt/2. Momenta written in the legends are spatial and in units of 2Tra/L. 
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FIG. 6: Effective value of Fy (At, At'; q^) at (Q, m) = (0, 0.015). 



with the point-to-all propagator where their sum At+At' is kept fixed. This is helpful to identify 
the plateau in Fy (At, At'; q^) as well as to stabilize the fit by increasing the number of available 
data. Fit results are summarized in Tables [Tl]-|Vl The statistical accuracy of Fy {q^ ) is typically 
3-5% at all of our simulated quark masses. Only the combination of our two largest momenta 
(|p|; Ip'I) = (V^) V^) has larger statistical error of about 10%. This is probably because precise 
determination of the ground state contribution is difficult due to very rapid damping of pion 
correlators with such large momenta. 
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TABLE II: Fit results for pion form factors at {Q, m) = (0, 0.050). 



IpI Ip'I |q| 




IpI Ip'I |q| 


„2 Fsiq'^) 

1 ^v[q ) ;.^(,.^^) 




1 1 
\/2 1 1 
V3 y/2 1 
\/2 \/2 


1 1.30(11) 
-0.1190(7) 0.7219(37) - 
-0.13763(68) 0.660(13) 0.918(37) 
-0.1436(23) 0.599(49) 0.907(90) 
-0.2083(18) 0.5908(45) 0.871(17) 


a/3 a/3 
V3 1 V2 
1 1 V2 
a/2 1 V3 
1 1 2 


-0.2865(95) 0.5132(83) 0.705(44) 
-0.2546(51) 0.522(17) 0.728(66) 
-0.3084(0) 0.4755(47) 0.759(21) 
-0.4461(7) 0.3726(72) 0.681(34) 
-0.6169(0) 0.3075(61) 0.605(33) 


TABLE III: Fit results for pion form factors at (Q, m) = (0, 0.035). 


IpI Ip'I |q| 




IpI Ip'I |q| 


„2 F,-(f,'^\ Fsil^) 
1 ^ Fsiqi,) 




1 1 
V2 1 1 
V3 \/2 1 
v'2 V2 


1 1.28(10) 

-0.1103(15) 0.7272(65) - 
-0.1350(17) 0.638(18) 0.902(43) 
-0.1417(27) 0.610(85) 0.890(87) 
-0.1873(44) 0.5964(75) 0.728(64) 


\/3 

V3 1 V2 
1 1 V2 
V2 1 V3 
1 1 2 


-0.2513(99) 0.5266(98) 0.673(69) 
-0.2459(54) 0.513(25) 0.783(66) 
-0.3084(0) 0.4518(80) 0.765(27) 
-0.4435(17) 0.358(11) 0.648(53) 
-0.6169(0) 0.2882(76) 0.546(31) 


TABLE IV: Fit results for pion form factors at {Q, m) = (0, 0.025). 


IpI Ip'I |q| 




IpI Ip'I |q| 


q ^v{q ) p^^gi^^ 




1 1 
V2 1 1 
Vs V2 1 
V2 V2 


1 1.29(13) 
-0.1010(27) 0.7327(88) - 
-0.1375(27) 0.627(21) 0.939(61) 
-0.1410(38) 0.550(60) 0.49(31) 
-0.1790(70) 0.600(10) 0.768(30) 


\/3 V3 
V3 1 V2 
1 1 V2 
V2 1 \/3 
1 1 2 


-0.237(15) 0.529(12) 0.650(35) 
-0.2489(83) 0.474(24) 0.56(11) 
-0.3084(0) 0.422(11) 0.686(42) 
-0.4460(27) 0.337(12) 0.548(50) 
-0.6169(0) 0.2561(87) 0.523(54) 


TABLE V: Fit results for pion form factors at {Q, m) = (0, 0.015). 


IpI Ip'I |q| 




IpI Ip'I |q| 


„2 p, („2\ Fs(q^) 




1 1 
\/2 1 1 
\/3 \/2 1 
\/2 \/2 


1 1.29(19) 
-0.0910(44) 0.727(13) - 
-0.1372(73) 0.645(34) 0.926(86) 
-0.1301(98) 0.629(74) 0.81(41) 
-0.162(21) 0.631(25) 0.805(44) 


\/3 VS 
V3 1 V2 
1 1 V2 
V2 1 \/3 
1 1 2 


-0.174(18) 0.574(14) 0.621(90) 
-0.227(11) 0.484(33) 0.80(12) 
-0.3084(0) 0.403(16) 0.643(95) 
-0.4456(73) 0.319(18) 0.599(81) 
-0.6169(0) 0.242(14) 0.65(13) 
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FIG. 7: Effective value of normalized scalar form factor Fs{At, At' ; q'^)/ Fs{At, At' ; q'^^^) at {Q,m) 
(0,0.050). 
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FIG. 8: Effective value of Fs{At,At';q^)/Fs{At,At';q^^i) at (g, m) = (0, 0.015). 
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B. Scalar form factor 



The scalar form factor normalized by the value at a certain momentum transfer q^^^ can be 
calculated from the following ratio similar to that for Fv{q'^) 

Fs{At,At';q^) ^ RsiAt,At';q^) 

FsiAt,At';ql,) Rs{At, At';ql,y ^ ' 

^ ^^IpUp'I fi.ed^,|p,, C..,,.,,(At;p)C..,,,,,(Ai';pO' ^ 



where 



C^f)(At,At';p,pO = C7^r)(At,At';p,p')-Citf)(At,At';p,p') (26) 

+ ci7;)(Ai,At';p,p') 
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FIG. 9: Normalized scalar form factor Fsiq^)/ Fs{q^^f) with and without the contributions of discon- 
nected diagram to Fsiq^)- We use the normalization value Fs{q^^f ) with the disconnected contributions 
in both data. 



is the three-point function with the flavor-singlet scalar operator. We note that C^^^^ — C^^^^ 
suffers from a severe cancellation: it is typically a subtraction of 0(1) quantities to extract 
their O(10~^) difference. Although this subtraction leads to a large statistical error in Fs{q'^), 
it is present only at |q| = 0, since C^-^^^ vanishes at nonzero |q|. In the following analysis, 
therefore, we use Fs{q'^) normalized at the smallest nonzero \q'^\ with |qref | = 1 rather than at 
5^ = 0. The effective value of the normalized scalar form factor Fs{/^t, At' ; q"^) / Fs{/^t, At' ; q'^^^) 
is plotted in Figs. [7] and El We summarize Fs{q^) / Fs{q^^^) determined from a constant fit to 
Fs{At,At';q^)/Fs{At,At';ql,) in Tables |TI]-|Vl 

Figure [9] compares Fs{q^) / Fs{q'^^^) to that without the contributions of the disconnected 
diagrams to Fs{q^). We observe a small but significant deviation between the two data, which 
implies the importance of the disconnected contributions. 

An accurate estimate of ^^(0) is useful for a precise determination of the scalar radius {r'^)s-, 
as it characterizes the q^ dependence of Fs{q^) near q^ = 0. Since the extraction of ^5(0) from 
the three-point function suffers from large statistical error due to the subtraction of the vacuum 
expectation value (VEV) of the scalar operator, we test an alternative calculation of ^5(0) 
through the Feynman-Hellmann theorem 

1 dM'^ 

Fsm = ~. (27) 

where the VEV subtraction is implicitly taken into account. Note that the overall factor 1/2 
is present in the RHS, since ra is the mass of two degenerate quark flavors. This and -FsC'Zref) 
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TABLE VI: Scalar form factor ^5(0) determined from Feymnan-Hellmann theorem Eq. ((27|) . and its 
value normalized by Fsiq^^.f) from ratio ([28| . The first error is statistical. The second is systematics due 
to the choice of the fit form for and uncertainty in Zs- 



m 


0.015 


0.025 


0.035 


0.050 


Fs{0) 


1.149(6)(22) 


1.162(18)(37) 


1.208(29)(47) 


1.319(44)(52) 


Fs{0)/Fs{qi,) 


1.413(70)(28) 


1.415(44)(44) 


1.338(38)(51) 


1.441(53)(58) 



determined from a ratio of pion correlators 

' AtAt'^00 ^ "^'^'^C^,,<^^^,(At;p)C,,,<^^^,(At + At';pO ^ ^ 

provide an alternative estimate of the normalized form factor Fs{0)/Fs(q^^^). We determine 
^5(0) from the chiral fit of and = 0.838(14)(3) presented in Ref. [sj. Results for Fs{0) 
and Fs (0) / Fs (q^^f) are summarized in Table IVll where the systematic error is estimated by 
changing the fitting form for and by taking account of the uncertainty in Z5. It turns 
out that, with the setup of our measurements, the Feynman-Hellmann theorem (j27p leads to a 
slightly smaller uncertainty in Fs{0)/ Fs{q'^^f) than that from the ratio (pl|) . We therefore adopt 
Fs{0)/Fs{q'^^{) in Table ED and Fsiq"^ j^O)/Fs{q'^^i) in Tables [TI]-|V] in the following analysis. 



C. Finite volume correction 



The finite volume effect could be significant at two smallest quark masses m = 0.015 and 
0.025, as the value of M.,^ L is less than 4. We estimate the finite volume correction (FVC) to 
the pion mass as presented in Ref. Qj. The FVC to the vector form factor Fy(g^) has been 
calculated within one- loop ChPT in Refs. 



451, 



4g] by replacing the loop integral by a discrete 



sum 



1 /"^ 

AFy(g') = ^ dx {h/2{L,xci,M^ + x{l - x)q^) - h/2{L,0,M^)} 



(29) 



where 



/,(L,Ak,M2 



(27r)3 {(k + Ak)2 + L3 ^ {(k + Ak)2 + M^Y 



-y 



(30) 
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FIG. 10: Three-point functions with vector current (left panels) and scalar operator (right panels) 
calculated in different topological sectors for m = 0.050. Top and bottom panels show data at the 
smallest and largest nonzero \q^\. 



This function can be evaluated through the elliptic 6 function 
Is{L,Ak,M^) 



47 



481] as 



(47r)3/2 



(31) 
(32) 



The FVC to the scalar form factor Fs{q^) is similarly evaluated as 



AFs(g2) 



B 

2^2 



dx 



2g2 + Ml 



/3/2(L,:Eq,M2 +x(l -x)g2) +/,/2(L,0,M2)| ^ (33) 



where we use i? = 2.10(45) GeV obtained from our analysis of the pion mass [5|. We find that 
/S.Fy{q^) has a mild cp' dependence and its magnitude is similar to or smaller than the statistical 
error: it is typically 3-5% at m = 0.015, and decreases down to < 1 % at m = 0.050. The FVC 
to the normalized scalar form factor Fs{q^) / Fg^q^^^ is below the statistical uncertainty even at 
our smallest quark mass due to a partial cancellation of FVCs in the ratio. We use Fy{q^) and 
Psi'f') / Fs{(l^ei) with the FVC included in the following analysis. 



D. Fixed topology effects 



The form factors listed in Tables [Tl]-|V] are extracted in the trivial topological sector, and are 
subject to effects of the fixed topology. The effects are known to be suppressed by the inverse 
space-time volume 1/V 



37 



and are systematically correctable [3S] . In order to confirm if the 
correction to our data is statistically insignificant, we repeat the calculation of Fy{q^) and Fs{q'^) 
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TABLE VII: Fit results for pion form factors at (Q, to) = (-2, 0.050). 



IpI 


Ip'I 


|q| 


9^ 






IpI 


Ip'I 


|q| 




Fv{q^) 




1 





1 


-0.1200(11) 


0.7145(46) 




V3 


1 


V2 


-0.2503(70) 


0.481(22) 


0.79(13) 


V2 


1 


1 


-0.1377(15) 


0.667(20) 


0.923(51) 


1 


1 


V2 


-0.3084(0) 


0.4640(65) 


0.738(29) 


73 


V2 


1 


-0.1416(33) 


0.662(56) 


0.32(34) 


V2 


1 


V3 


-0.4461(15) 


0.365(11) 


0.694(46) 


V2 





n/2 


-0.2102(40) 


0.5910(75) 


0.829(23) 


1 


1 


2 


-0.6169(0) 


0.2962(89) 


0.635(55) 


V3 





V3 


-0.281(12) 


0.5113(97) 


0.785(65) 















TABLE VIII: Fit results for pion form factors at (Q, to) = (-4, 0.050). 



IpI 


Ip'I 


|q| 


q' 


Fv{q^) 


Fsiq^) 
Fs(ql,} 


IpI 


Ip'I 


|q| 


q' 


Fv{q^) 


Fs{q^) 
Fs(ql,) 


1 





1 


-0.1204(9) 


0.7234(58) 




^/3 


1 


V2 


-0.2660(42) 


0.561(22) 


0.94(11) 


V2 


1 


1 


-0.1357(11) 


0.641(21) 


0.899(49) 


1 


1 


V2 


-0.3084(0) 


0.4810(83) 


0.846(29) 


V3 


V2 


1 


-0.1493(16) 


0.699(76) 


0.39(32) 


V2 


1 


^/3 


-0.4442(11) 


0.3714(96) 


0.754(57) 


V2 





V2 


-0.2062(29) 


0.5875(48) 


0.872(26) 


1 


1 


2 


-0.6169(0) 


0.3143(86) 


0.603(42) 


V3 





x/3 


-0.3107(78) 


0.5392(87) 


0.820(52) 















in non-trivial topological sectors with Q = —2 and —4 at m = 0.050. Note that, LO ChPT predicts 
(Q2) = mSy/2~8 with our estimate of the chiral condensate S^(2 GeV) = (0.236(1:^^) GeV)^ 



[5|]. It is therefore not necessary to simulate topological sectors with |Q| » \/{Q^ ~ 3 at this 
quark mass. 

We compare three-point functions calculated in the different topological sectors in Fig. [TUl 
where no systematic deviation among the data is observed. This is also the case for Fy^siQ"^) 
summarized in Tables IVTIl and rVIIII The form factors at Q = 0, —2 and —4 are consistent with 
each other within two standard deviations as shown in Fig. [TTJ Although the effect of the fixed 
topology is likely below our statistical accuracy, we take the spread in Fv,s{<f ) as a systematic 
error at m = 0.050 and an uncertainty of the same magnitude is assumed at m < 0.050 in the 
following analysis. 
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FIG. 11: Vector (left panels) and scalar form factors (right panels) calculated in different topological 
sectors for m = 0.050. 



IV. PARAMETRIZATION OF THE DEPENDENCE 
A. Vector form factor 

The vector form factor Fy{q'^) is plotted as a function of q'^ in Fig. [T2j Its dependence 
turns out to be close to the expectation from the vector meson dominance (VMD) hypothesis 



where Mp is the vector meson mass calculated at the same quark mass. Then, we assume that 
the small deviation due to higher poles or cuts can be well parametrized by a polynomial form. 
We therefore fit our data to the following form 

^v{q'^) = ^TTTo + «v^,i 9^ + ay,2 (g^)^ + ay,3 (g^)^ (35) 



1 - g2/M 



in order to extract the charge radius {r )v and the curvature cy 



{r^)v = 6 



dFviq') 



d{q^) 



cv 



d'Fviq') 



(36) 

g2=0 



a(g2)2 

The fit curve is plotted in Fig. [12] and numerical results are summarized in Table llXi The fit 
describes our data reasonably well, and results for {r'^)v and cy do not change significantly by 
including or excluding the cubic term av^siq"^)^- In the following analysis, we employ results 
from the parametrization with the cubic term. 

One of the main purposes of this work is to investigate whether the q^ dependence of our data 
can be described by two- loop ChPT [isl. [l^. Figure [T3l shows contributions to Fv{q^) from each 
order (q'^)" of a Taylor expansion of Eq. (|35p . We find that 0{q^) and higher order contributions 
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q [GeV^] q [GeV"] 

FIG. 12: Vector form factor Fv^q^) as a function of q^. Solid and dotted lines show the fit curve ([55)1 
and its error. The dependence from the vector meson dominance model is shown by the dashed line. 



TABLE IX: Parametrization Eq. ([55]) for vector form factor Fv{q^)- Results for the vector radius {r'^)v 
and curvature cy in lattice units are also listed. 



m 


X^/d.o.f. 


av,i 


ay,2 


a.v,3 


{r')v 




0.050 


1.8(0.8) 


0.187(33) 


0.181(55) 




18.72(20) 


8.786(55) 


0.050 


1.7(0.9) 


0.116(64) 


-0.22(26) 


-0.48(28) 


18.30(38) 


8.38(26) 


0.035 


1.3(0.5) 


0.136(60) 


0.124(96) 




20.15(36) 


10.51(10) 


0.035 


1.2(0.7) 


0.02(11) 


-0.50(44) 


-0.73(46) 


19.48(66) 


9.88(44) 


0.025 


1.9(0.6) 


0.034(90) 


-0.03(14) 




21.69(54) 


12.78(14) 


0.025 


1.7(0.7) 


-0.14(16) 


-0.98(59) 


-1.10(56) 


20.64(98) 


11.83(59) 


0.015 


0.8(0.8) 


-0.09(12) 


-0.15(19) 




22.51(72) 


14.58(19) 


0.015 


0.5(0.8) 


-0.35(25) 


-1.56(96) 


-1.60(95) 


20.9(1.5) 


13.18(96) 



to Fviq"^) are sufficiently small only below ~0.3 GeV^, which is however around our smallest 
value of We therefore do not use the parametrization based on ChPT in this study. Note 



20 



-4.0 







m = 0.025 






— total 






— Oiq) 






— 0{q) 


- 




' ' O(q') + higher - 










q- [GeVl 



FIG. 13: Contributions in the expansion of Fv{q^)—1 at r7i = 0.025. Thin solid, dashed and dot-dashed 
hnes show 0{q^), 0{q*) and higher contributions, whereas the thick sohd hne is their totaL 



TABLE X: Single pole fit ((37)) . We also list the p meson mass at simulated quark masses. 



m 


xVd.o.f. 






cv 


Mp 


0.050 


4.3(1.1) 


0.5431(50) 


20.34(38) 


11.49(43) 


0.5839(28) 


0.035 


2.3(0.7) 


0.5270(80) 


21.61(65) 


12.97(78) 


0.5570(31) 


0.025 


2.9(0.8) 


0.511(11) 


23.0(1.0) 


14.6(1.3) 


0.5285(43) 


0.015 


0.7(0.7) 


0.520(16) 


22.2(1.4) 


13.7(1.7) 


0.5104(55) 



that such large higher order contributions are unavoidable unless | | <C Mp , because the VMD 
form is a good approximation of Fy(q'^). 

We also test a single pole ansatz often used in the previous studies 

^ T^Ti^; '''' 

As summarized in Table 1X1 this fit tends to give a slightly higher ami (r^)y than those from 
Eq. (|35]) . This may suggest that it is difficult to describe our precise data of FY{q'^) in the whole 
region of g^[GeV^] € [—1.7,0] by a simple pole-dominance form. 



B. Scalar form factor 



Due to the lack of knowledge about the scalar resonances at the simulated quark masses, we 
test a generic polynomial form up to the quartic order 

Fs{q^) = Fs{0) { 1 + as,i + as,2 {q^f + as,3 (q^f + as,4 {qY} (38) 
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FIG. 14: Normalized scalar form factor Fs (q^) / Fs (q^^f) as a function of q^. Solid and dotted lines show 
the cubic fit and its error. 



TABLE XI: Parametrization Eq. p8|) for scalar form factor Fs{q'^)- Results for the vector radius {r'^)s 
and curvature cs in lattice units are also listed. 



m 


xVd.o.i. 


as,i 


as,2 


as,3 


as,4 


{r')s 


Cs 


0.050 
0.050 


1.3(1.0) 
1.3(1.1) 


3.04(31) 
3.51(49) 


6.6(1.3) 
10.5(3.5) 


5.1(1.4) 
1.6(9.2) 


9.6(7.7) 


18.2(1.9) 
21.1(3.0) 


6.6(1.3) 
10.5(3.5) 


0.035 
0.035 


1.8(1.0) 
1.9(1.2) 


2.79(31) 
2.60(51) 


5.6(1.4) 
3.5(4.8) 


4.3(1.6) 

-3(15) 


-6(13) 


16.8(1.9) 
15.6(3.1) 


5.6(1.4) 
3.5(4.8) 


0.025 
0.025 


1.9(1.1) 
1.9(1.3) 


3.37(32) 
2.97(53) 


6.1(1.6) 
1.8(4.6) 


3.6(1.8) 

-9(13) 


-12(11) 


20.2(1.9) 
17.8(3.2) 


6.1(1.6) 
1.8(4.6) 


0.015 
0.015 


1.5(1.0) 
1.7(1.1) 


3.51(51) 
3.0(1.0) 


6.6(2.7) 
1.0(9.1) 


3.6(3.3) 

-14(26) 


-16(22) 


21.0(3.1) 
18.1(6.0) 


6.6(2.7) 
1.0(9.1) 



to parametrize the dependence of the scalar form factor Fs{q'^). We observe that the cubic 
('^S',4 = 0) fit also describes our data reasonably well as seen in Fig. Fit results summarized 
in Table IXII show that the inclusion of the quartic correction does not change the value of 
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and the result for the scalar radius {r'^)s = Qo,s,i significantly. However, such a stability against 
the choice of the parametrization is much less clear in the curvature cs = 05,2 due to its large 
uncertainty. From these observations, we only use results for {r^)s in the following analysis, and 
leave a precise determination of C5 for future studies. 



CHIRAL EXTRAPOLATION 



A. Fit based on one-loop ChPT 



Since the form factors Fy^siff') are independent of at LO in ChPT, the chiral expansion 
of the radii {r^)v,s starts from the one-loop order of ChPT. We first compare our lattice results 



with the one-loop expressions 
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(r2)y = L^(l+6iV/n L^ln 

NF^ NF^ 



(39) 



Ml 



(40) 



where A^ = (47r)^, and F is the decay constant in the chiral limit. We adopt the normalization 
of the decay constant Ft^ = 92 MeV at the physical quark mass. The renormalization scale fj, is 
set to AttF in our analysis. At this order of the chiral expansion, F is the only LEC appearing 
in the M^^ dependent terms. We fix this important parameter to F = 79.0(lt2'6) MeV, which 
has been determined from our detailed analysis of the pion mass and decay constant [5]. Each 
one of Eqs. (j39p and (j40p . therefore, has a single fit parameter, namely LECs /g or in their 
constant term. 

We find that the NLO fits of lattice data are not quite successful as seen in Fig. [15] and 
Table Kni While the data of {r^)v can be fitted with reasonable x'^/d-o.f. 0.14, the value 
extrapolated to the physical quark mass 

{r^)y = 0.3637(43) fm^ is significantly smaller than 
the experimental value 0.437(16) fm^ based on Nf = 2 ChPT [16] and 0.452(11) fm^ quoted by 
PDG [49]. As for the scalar radius, the one- loop formula fails to reproduce our data of {r'^)s as 
indicated by the quite large value of x^/d.o.f. ~9 : the data have a mild quark mass dependence 
in contrast to the 6 times enhanced chiral logarithm compared to {r'^)v- 

This failure of the NLO fits is not due to our choice of F. If F is treated as a free parameter, 
the fit to {r'^)s results in an unacceptably large value F ~ 200 MeV to achieve reasonable 
X^/d.o.f. < 1, whereas our data of {r'^)v favor F~80 MeV. Thus we can not make a consistent 
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FIG. 15: Chiral fit of (r^)y (left panel) and {r'^)s (right panel) using one-loop ChPT formulae. Filled 
squares are the lattice data and the value extrapolated to the physical point. In the left panel, we also 
plot the experimental value {r'^)v — 0.437(16) fm^ from an analysis based on Nf — 2 ChPT 16] (open 
circle) and 0.452(11) fm^ quoted by Particle Data Group |49| (star). The star symbol in the right panel 
represents (r^)5 = 0.61(4) fm^ obtained from an indirect determination through tttt scattering (5o| . 



TABLE XII: Results of chiral extrapolations of (r^)y and {r'^)s using one-loop ChPT formulae Eqs. ([39]) 
and (|40p . The radii extrapolated to the physical point are also listed. 



xVd-oi. 


Fe X 10^ 




xVd.o.f. 


11 X 10^ 




0.14 


-6.59(12) 


0.3637(43) 


9.0 


2.94(39) 


0.797(15) 



analysis. We note that we have experienced a similar situation, namely (r^)y smaller than 
experiment and a small quark mass dependence of {r'^)si in our previous study with a different 
lattice action [23.], though simulated quark masses are heavier than those in this study. 

We also note that the lattice data for the curvature cy largely deviate from its NLO ChPT 
expression as shown below. From all of these observations we conclude that the chiral behavior 
of the pion form factors in the quark mass region ms/Q-ms/2 are not well described by NLO 
ChPT. 



B. Fit based on two-loop ChPT 

In Ref. [sj], we observed that NNLO contributions are important to reliably extract LECs 
from the pion mass and decay constant in our simulated region of the quark mass. Therefore 
there exists a possibility that NNLO contributions become also important to describe the chiral 
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TABLE XIII: Phcnomenological estimates of LECs used in this paper. We note that O(p^) couphngs 
^^vs\ Wc\ based on a resonance saturation hypothesis. 



Ref.[51] 


Ref.[50] 


Ref.fW] 


Ref.[16] 


F[MeY] 


h h 


^3 


k r'v^r X 10"* r\,c X 10^ x 10^ 


86.2(5) 


-0.36(59) 4.31(11) 4.39(22) 


2.9(2.4) 


16.0(9) -2.5 2.6 -0.3 



behaviour of the radii. The two- loop expressions of {r'^)v,s and cy are given by 
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where we use a hnear combination I12 — ^2/^ instead of l\, since the former is convenient 
for our chiral extrapolation (see below). The analytic terms containing '^{V5}{rc} ^'spresent 
contributions of tree diagrams with vertices from the 0(p^) chiral Lagrangian. 

Before fitting lattice data to these expressions, one can get some idea about the significance 
of the NNLO contributions by using phenomenological estimates of LECs. A collection of recent 
estimates is shown in Table IXlIll where the LECs in the 0{p^) chiral Lagrangian are denoted 
by the scale invariant convention /j defined by 



2N 
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li + In 
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Figure [T6l shows the expected Ml dependence of (r^)y and {r'^)s from Eqs. ()4ip and (j42p with 
the phenomenological estimates of LECs. The individual contributions from NLO and NNLO are 
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FIG. 16: Chiral behavior of radii based on two-loop ChPT with phenomenological estimates of LECs 
hsted in Table IXIIII The dashed and dot-dashed lines represent contributions at NLO and NNLO, 
whereas the solid line is their total. We note that r^,, ei r i from the resonance saturation are taken as 

{V,S}.,{i\c} 

the renormalized LECs at the resonance mass scale in this plot. 



also plotted. This analysis suggests that the NNLO contributions could significantly modify the 
chiral behavior of the radii in our simulated quark masses. However, we note that ^^ygy c} 
poorly known and these in Table [XlIIl are determined from a resonance saturation hypothesis. A 
chiral extrapolation of lattice data with the two-loop formulae is therefore important to confirm 
the significance of the NNLO contributions and to resolve the failure of the one-loop fit. 

The curvature cy characterizes the 0{q'^) dependence of Fy{q^), and therefore requires the 
NNLO terms to describe its logarithmic dependence on the quark mass as well as a constant 
term. Near the chiral limit it has a divergent term of the form ~ 1/(F^M^), which comes from 
non-analytic NLO contributions in Fy[q^). Since the divergent term 1/(F^M^) is significant 
only below the physical pion mass, the analysis of the lattice data for cy requires the NNLO 
contributions. 



We extend our analysis to two-loop ChPT as already outlined in our previous report 



The curvature cy is included into our chiral extrapolation to obtain an additional constraint 
on LECs. Both of {r'^)y and cy depend on l\ and I2 only through the linear combination ?^2) 
and the complicated two-loop expressions for (r^)y and cy involve only four free parameters /g, 
^12' "^Vr ^'^'^ "^Vc by choosing M.^/(47rF)^ as an expansion parameter. Therefore we first try 
a simultaneous fit to {r'^)y and cy to check that the chiral behavior of our data is described 
by two-loop ChPT. Fit curves are plotted in Fig. [T7] and numerical results are summarized in 
Table IXTVl The fit leads to an acceptable value of x^/d.o.f. ~ 0.7, and the relevant LECs are 
determined with a reasonable accuracy. We note that this fit is based only on ChPT without 
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TABLE XIV: Results of simultaneous chiral fit to (r^)y and cy based on two-loop formulae Eqs. (|¥T|) 
and (gni). 





11 X 10^ 


X 10^ 


r\,,. X 10^ 


rVc X 10'' 




cy[GeV~'*] 


0.70 


-8.48(87) 


-3.3(1.1) 


-0.77(69) 


3.97(20) 


0.411(26) 


3.26(21) 



additional assumptions. The extrapolated values of (r^)y and cy are consistent with recent 
phenomenological determinations using experimental data of Fv{q^) id. l52l. 



In Sec. IIVI we observe that the 0{q^) contribution to Fv{q^) is not small in the simulated 
region of the momentum transfer > 0.3 GeV^ (g^/(47rF)^ > 0.3). The dependence of our 
data is therefore parametrized by the generic polynomial forms (j35p and (j38p instead of those 
based on ChPT. This is not in contradiction with the successful chiral extrapolation of 
and Cy: since we explore small pion masses < 0.3 GeV^ (M^/(47rF)^ < 0.3), the quark mass 
dependence of (r^)y and cy is described by the two-loop ChPT formulae. 

The inclusion of {r'^)s into the simultaneous chiral fit introduces additional four free param- 
eters ^4 'iiid r^^, and we need to fix some of them to obtain a stable fit. Since ^2 
/g appear only in the NNLO terms and have been determined with a reasonable accuracy from 
phenomenology or lattice studies, we use a phenomenological estimate = 4.31(11) js^ and a 
lattice estimate /3 = 3.38(56) from our analysis of the pion spectroscopy Q]. We treat r^^, and 
l\ as free parameters because of poor knowledge on the former and in order to examine the 
consistency of the latter with that determined from F-,^. The fit curves are shown in Fig. [18] 
and numerical results are summarized in Table IXVl We observe that i) the results in the vector 
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TABLE XV: Results of simultaneous chiral fit to {r'^)v.s and cy based on two-loop formulae Eqs. (|4T|) 
(ga) and (gSl). 



xVd.o.f. ;g X f 0^ II X 10^ 2 X 10^ rl^r x 10^ r^>_^ x 10^ rj^ x lO'^ 



{r^)v[fm'] cv[GeV-4] 



.[fm^ 



0.68 -8.41(76) 1.1(3.2) -3.10(90) -1.0(1.1) 4.00(17) 1.74(36) 



0.409(23) 3.22(17) 0.617(79) 



channel, namely cy and the relevant LECs, do not change significantly by including {r'^)s 

into the chiral fit, and ii) this fit describes the mild quark mass dependence of {r'^)s with an ac- 
ceptable value of x^/d.o.f. ~0.7. From Fig. [TSl we find that the net NNLO contribution to {r'^)s 
is larger than NLO around our largest quark mass. This is due to an accidental cancellation 
between the constant and logarithmic terms at NLO. 

We take results from this simultaneous fit as our best estimate of {r'^)v,Sj cy. The systematic 
error due to the chiral extrapolation is estimated by excluding the data at the largest quark 
mass from the fit. The shifts in the radii and curvature turn out to be at the 1 a level: namely, 
the extrapolation is stable against variation of the fit range. In order to estimate systematics 
due to the fixed topology and the use of the Feynman-Hellmann theorem to estimate Fs{0), we 
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TABLE XVI: Statistical and systematic errors of radii and curvature. We list the magnitude of the 



errors relative to the central values. 





Stat. 


sys. (total) 


chiral fit 


input(LECs) 


input (ro) 


fixed Q 


Feynman-Hellmann 


finite a 




6% 


9% 


5.3% 


0.2% 


6.1% 


2.2% 


1.4% 


3.0% 




13% 


11% 


4.0% 


2.2% 


7.7% 


2.4% 


4.6% 


3.0% 


cv 


5% 


11% 


2.7% 


0.4% 


9.3% 


2.3% 


3.8% 


3.0% 



repeat the whole analysis by using Fy^s{<f) shifted by their systematic uncertainties discussed in 
Sees. HITB] and UlI D[ Systematic uncertainties due to the choice of the inputs, namely LECs F, 
I2, ^3 for the chiral extrapolation and tq = 0.49 fm to fix the scale, are estimated by shifting each 



55|, 



561 ] ■ In addition, 



LEG by its uncertainty and by using a recent lattice estimate rg = 0.47 fm 
a discretization error estimated by a naive order counting 0((aAQCD)^) ~ 3% is also taken into 
account. The magnitude of these systematic uncertainties are summarized in Table IXVTl Our 
final results for the radii and curvature at the physical point are 

{r^)v = 0.409(23)(37) fm^, (46) 
{r^)s = 0.617(79)(66) fm^, (47) 
cv = 3.22(17)(36) GeV"^, (48) 

where the first error is statistical and the second represents systematic errors added in quadra- 
ture. These results are in good agreement with their phenomenological values. The largest 
systematic error arises from the use of the input rg to fix the scale. This can be removed by a 
better choice, such as f-j^, in future studies in three-flavor QCD. 

We estimate the relevant LECs in a similar way to that for {r'^)v,s cy, and obtain 

k = 11.9(0.7)(1.0), (49) 

k = 4.09(50) (52), (50) 

h-h = -2.9(0.9)(1.3), (51) 

r;>_^ = -1.0(1.0)(2.5) X 10~^ (52) 

r|>^c = 4.00(17)(64) x 10"^ (53) 

r^^^ = 1.74(36)(78) x 10"^ (54) 

Our estimate of /g is slightly smaller than those from two-loop ChPT analyses of experimental 
data; ^6 = 16-0(0.9) from Fy [16] and 15.2(0.4) from r and vr decays [57|]. This is partly due to 
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the deviation of F between our lattice determination |5l and two-loop ChPT 1511. We note that 
I4, is consistent with our lattice determination /4 = 4.12(56) from F.,^ 5] and a phenomenological 
estimate 4.39(22) jsol. Here we list I1 — I2 instead of 1^2-, which is scale invariant and equals to 
(/i — /2)/(6A^) in the convention of Eq. (j44l) . Although I1—I2 has a large uncertainty of ~ 50%, 
it is consistent with —4.67(60) from phenomenology. Note that we set the renormalization scale 
^ = At:F, and our fit favors small values of the order of 10~^-10~^ for the renormalized 0{p^) 
couplings jrc} ^^^^ scale. Finally, we emphasize that the simultaneous fit only to (r^)y 

and cy without any phenomenological input leads to consistent results for the vector channel, 



namely for (r^)y, cy, 1^, h — h and r 



V,{r,c}- 



VI. CONCLUSIONS 



In this article, we present a lattice calculation of the pion form factors Fy^sil'^) in two- flavor 
QCD. We study the chiral behavior of the radii {r'^)v,s and the curvature cy based on ChPT up 
to two loops. We employ the overlap quark action, which has exact chiral symmetry and thus 
provides the cleanest framework for the study of chiral behavior. Otherwise, the explicit chiral 
symmetry breaking of conventional lattice actions could make studies based on two-loop ChPT 
substantially complicated. Another salient feature of this work is that, for the first time, Fs{q'^) 
is evaluated including the contributions of the disconnected diagrams through the all-to-all quark 
propagators. 

Our detailed analysis based on ChPT reveals that two-loop contributions are important 

to describe the chiral behavior of {r'^)v.S and cy at our simulated quark masses, which are 

comparable to those in recent unquenched simulations. Through the chiral extrapolation at two 

loops, we obtain {r'^)v,s and cy at the physical point, which are consistent with experiment. 

We also obtain estimates of O(p^) and 0{p^) LECs, and confirm that Fs(q^) and F^^ lead to 

consistent results for [4 as suggested long ago jl4| . 

I — I 

As we already outlined in Ref . [28] , the curvature cy is useful to stabilize the two- loop chiral 
fit. The single pole ansatz Eq. (j37p . which has been commonly used in previous studies, may not 
be suitable to estimate cy, since it simply assumes a relation cy = ((r^)y/6)^ at simulated quark 
masses. In this work, instead, we employ a generic form up to cubic corrections to parametrize 
the q'^ dependence of Fy^siq^) thanks to the precise determination of Fy^siq^) through the all- 
to-all propagator. Dispersive analyses of the q^ dependence, model independent information 
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of scalar resonances, and the twisted boundary condition [58] already used in Refs.Q, [^^J are 
interesting subject and technique for a better control of this parametrization in future studies. 

Extending this study to three-flavor QCD is important for more realistic comparison with 
experiment. We have already started simulations in three-flavor QCD with the overlap action 
; measurements of pion correlators using all-to-all propagators are in progress. 

Finally, it is expected from our chiral extrapolation that {r'^)s shows a strong quark mass 
dependence due to the one-loop chiral logarithm below M-,^ ~ 250 MeV. Pushing simulations 
toward such small quark masses is an interesting subject in future studies for a direct observation 
of the one-loop logarithm, although it is very challenging. 
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